##########################################################################################################################################################
### Setup

### Source setup
source("code/1_modeling_setup.R")



########################################################################################################################################################################
### Table 1: Divided government summary statistics

### Load and format divided government data

# Load, add divided indicator and reformat state name
pty <- assignLoad("data/partyControl.RData") %>% 
  mutate(divided = ifelse(state_control=='Divided', 1, 0),
         state = state.abb[match(state, state.name)])

# Format indicator for whether there was a change to divided government
div_change <- pty %>%
  dplyr::select(state, year, divided) %>% 
  mutate(divided_lag = lag(divided, order_by=state, default=NA)) %>%
  mutate(change_to_div = ifelse(divided==1 & divided_lag==0, 1, 0),
         change_from_div = ifelse(divided==0 & divided_lag==1, 1, 0)) %>%
  mutate(divided_lag = ifelse(is.na(divided_lag), 0, divided_lag),
         change_to_div = ifelse(is.na(change_to_div), 0, change_to_div))
div_change_merge <- div_change %>%
  dplyr::select(state, year, change_to_div, change_from_div)

# Indicator for whether a state ever had divided government
div_any <- pty %>%
  group_by(state) %>%
  dplyr::summarise(
    div_any = ifelse(any(divided==1), 1, 0)
  ) %>%
  ungroup()

# Join div_lab and div_any back to pty
pty <- pty %>%
  left_join(., div_change_merge, by=c("state", "year")) %>%
  left_join(., div_any, by="state")

# Fix one weird case (AK 2009) manually
pty$change_to_div[is.na(pty$change_to_div)] <- 0


### Divided government summary table

# 1) Number and proportion Republican-held chambers by year
# 2) Number and proportion of divided governments
# 3) Number and proportion of states changing to divided government
# 4) Number and proportion of states changing from divided government
div_summ <- pty %>%
  mutate(rep_hs_maj = ifelse(hs_maj_party=='D' | is.na(hs_maj_party), 0, 1),
         rep_sen_maj = ifelse(sen_maj_party=='D' | is.na(sen_maj_party), 0, 1)) %>%
  group_by(year) %>%
  dplyr::summarize(
    # Number of chambers
    num_states = n(),
    # Number and proportion of chambers held by Republicans
    num_rep_cham = sum(rep_hs_maj) + sum(rep_sen_maj),
    prop_rep_cham = num_rep_cham/(num_states*2)*100,
    # Number and proportion of states with divided government
    num_div = sum(divided),
    prop_div = num_div/num_states*100,
    # Number and proportion of states changing to divided government
    num_to_div = sum(change_to_div),
    prop_to_div = num_to_div/num_states*100,
    # Number and proportion of states changing from divided government
    num_from_div = sum(change_from_div),
    prop_from_div = num_from_div/num_states*100
  ) %>%
  ungroup() %>%
  round()

# Prettify for putting into a table
div_summ_tab <- div_summ %>%
  filter(year %in% 2010:2014) %>%
  dplyr::rename(Year=year) %>%
  mutate(`Republican Control` = paste0(num_rep_cham, ' (', prop_rep_cham, '%)'),
         `Divided Government` = paste0(num_div, ' (', prop_div, '%)'),
         `Change to Divided` = paste0(num_to_div, ' (', prop_to_div, '%)'),
         `Change to Unified` = paste0(num_from_div, ' (', prop_from_div, '%)'),
         Year = char(Year)) %>%
  dplyr::select(Year, `Republican Control`, `Divided Government`, `Change to Divided`, `Change to Unified`)

# Save
div_tab <- xtable(div_summ_tab)
print(div_tab, file='results/divided_government_by_year.tex')


### Add relevant variables to bill and votes

# Subset
div_merge <- dplyr::select(pty, c(state, year, div_any, change_to_div, change_from_div))

# Format change_to_div for plotting 
div_merge <- div_merge %>%
  mutate(change_to_div_plot = ifelse(change_to_div==1, 'Changed to Divided Govt', 'No Change') %>% as.factor() %>% relevel(ref='No Change'))

# Merge
bill <- left_join(bill, div_merge, by=c('state', 'year'))
votes <- left_join(votes, div_merge, by=c('state', 'year'))


### Make indicator for state-years that changed to divided government *and* the year before

# Get state years
todiv_styr <- filter(div_change, divided==1 & divided_lag==0 & year %in% 2010:2014) %>%
  mutate(state_year = paste0(state, '_', year)) %>%
  pull(state_year)

# Add the year before
prev_todiv_styr <- c()
for(styr in todiv_styr){
  prev_yr <- str_extract(styr, '[0-9]{4}') %>% num() %>% `-`(1) %>% char()
  result <- str_replace(styr, '[0-9]{4}', prev_yr)
  prev_todiv_styr <- c(prev_todiv_styr, result)
}

# Merge
div_styrs <- c(todiv_styr, prev_todiv_styr) %>%
  str_delete("^(NY|IA|ME|NM|NJ)") # Remove observations where we don't have coverage pre-post for each state
# Votes data is also missing CO, NH, OR



##########################################################################################################################################################
### Figure 1: Correlation matrix for all variables

### Correlation matrix

# Get variable names
votes %>%
  dplyr::select()

# Format ivMat
ivMat <- votes %>%
  dplyr::select(veto_override_threshold, lineitemveto, yea_prop, min_yea_prop, unity_vote, 
                seat_ratio, term_diff, divided3, committee_sponsor,
                maj_sponsor, num_cs, maj_share, year, chamber) %>%
  mutate(divided3 = ifelse(divided3=='Divided', 1, 0),
         house = ifelse(chamber=='H', 1, 0),
         maj_sponsor = car::recode(maj_sponsor, "'Majority'=1; 'Minority'=0") %>% num(),
         unity_vote = num(unity_vote),
         year = num(year),
         lineitemveto = num(lineitemveto),
         committee_sponsor = num(committee_sponsor)) %>%
  dplyr::select(-chamber)

# Create correlation matrix
corMat <- cor(ivMat, use='pairwise.complete.obs')

# Format names for plotting
rownames(corMat) <- colnames(corMat) <- c('Veto Override Threshold', 'Line Item Veto', 'Proportion Voting Yea', 'Proportion Minority Voting Yea',
                                          'Party Unity Vote', 'Seat Ratio', 'Difference in Term Length', 
                                          'Divided Government', 'Committee Sponsor', 'Majority Sponsor', 'Number of Cosponsors', 
                                          'Majority Seat Share', 'Year', 'House')

# Plot correlation matrix, save
pdf(file='~/Dropbox/Projects/Agenda Access in States/Results/model_variables_correlations.pdf', width=10, height=10)
corrplot(corMat, method='circle', tl.col='black', type='upper')
dev.off()



###############################################################################################################################################################
### Main regression tables (Tables 3-6)

### Models

# Interchamber IVs, pr(vote)
mods_vote_divided_inter_int <- vector(mode='list', length=2)
mods_vote_divided_inter_int[[1]] <- glm(voted_sles ~ veto_override_threshold*divided3 + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))
mods_vote_divided_inter_int[[2]] <- glm(voted_sles ~ lineitemveto_plot*divided3 + veto_override_threshold + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))

# Interchamber IVs, pr(law)
mods_law_divided_inter_int <- vector(mode='list', length=2)
mods_law_divided_inter_int[[1]] <- glm(law_sles ~ veto_override_threshold*divided3 + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))
mods_law_divided_inter_int[[2]] <- glm(law_sles ~ lineitemveto_plot*divided3 + veto_override_threshold + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))

# Interchamber IVs, prop(yea)
mods_yea_divided_inter_int <- vector(mode='list', length=2)
mods_yea_divided_inter_int[[1]] <- lm(yea_prop ~ veto_override_threshold*divided3 + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))
mods_yea_divided_inter_int[[2]] <- lm(yea_prop ~ lineitemveto_plot*divided3 + veto_override_threshold + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))

# Interchamber IVs, prop(yea|min)
mods_yea_min_divided_inter_int <- vector(mode='list', length=2)
mods_yea_min_divided_inter_int[[1]] <- lm(min_yea_prop ~ veto_override_threshold*divided3 + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))
mods_yea_min_divided_inter_int[[2]] <- lm(min_yea_prop ~ lineitemveto_plot*divided3 + veto_override_threshold + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))

# Interchamber IVs, pr(unity)
mods_unity_divided_inter_int <- vector(mode='list', length=2)
mods_unity_divided_inter_int[[1]] <- glm(unity_vote ~ veto_override_threshold*divided3 + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness))
mods_unity_divided_inter_int[[2]] <- glm(unity_vote ~ lineitemveto_plot*divided3 + veto_override_threshold + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness))

# Bicameralism IVs, pr(vote)
mods_vote_divided_bicam_int <- vector(mode='list', length=2)
mods_vote_divided_bicam_int[[1]] <- glm(voted_sles ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))
mods_vote_divided_bicam_int[[2]] <- glm(voted_sles ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))

# Bicameralism IVs, pr(law)
mods_law_divided_bicam_int <- vector(mode='list', length=2)
mods_law_divided_bicam_int[[1]] <- glm(law_sles ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))
mods_law_divided_bicam_int[[2]] <- glm(law_sles ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))

# Bicameralism IVs, prop(yea)
mods_yea_divided_bicam_int <- vector(mode='list', length=2)
mods_yea_divided_bicam_int[[1]] <- lm(yea_prop ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))
mods_yea_divided_bicam_int[[2]] <- lm(yea_prop ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))

# Bicameralism IVs, prop(yea|min)
mods_yea_min_divided_bicam_int <- vector(mode='list', length=2)
mods_yea_min_divided_bicam_int[[1]] <- lm(min_yea_prop ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))
mods_yea_min_divided_bicam_int[[2]] <- lm(min_yea_prop ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))

# Bicameralism IVs, pr(unity)
mods_unity_divided_bicam_int <- vector(mode='list', length=2)
mods_unity_divided_bicam_int[[1]] <- glm(unity_vote ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness))
mods_unity_divided_bicam_int[[2]] <- glm(unity_vote ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness))


### Tables

# Line-item veto
line_mods <- list(mods_vote_divided_inter_int[[2]], mods_law_divided_inter_int[[2]], mods_yea_divided_inter_int[[2]], 
                  mods_yea_min_divided_inter_int[[2]], mods_unity_divided_inter_int[[2]])
stargazer(line_mods, keep.stat = 'n', font.size='tiny', column.sep.width='1pt', digits=2, 
          dep.var.caption='', header=FALSE, label='t.line.inter', model.names=FALSE, model.numbers = FALSE,
          dep.var.labels=c("Reached Vote", "Became Law", "Prop Yeas", "Prop Minority Yeas", "Unity vote"),
          title='Estimates and standard errors for models showing the relationship between line-item veto and our dependent variables.',
          covariate.labels = c("Line-Item Veto", "Divided Government", "Veto Override Threshold", "Majority Sponsor", 
                               "Number of Cosponsors", "Majority Seat Share", "Party Difference", "Legislative Professionalism",
                               "Committee Sponsor", "Majority Unity", "Governor's Party can Override",
                               2011:2014, "Line-Item Veto * Divided Government"))

# Veto override threshold
veto_mods <- list(mods_vote_divided_inter_int[[1]], mods_law_divided_inter_int[[1]], mods_yea_divided_inter_int[[1]], 
                  mods_yea_min_divided_inter_int[[1]], mods_unity_divided_inter_int[[1]])
stargazer(veto_mods, keep.stat = 'n', font.size='tiny', column.sep.width='1pt', digits=2, 
          dep.var.caption='', header=FALSE, label='t.veto.inter', model.names=FALSE, model.numbers = FALSE,
          dep.var.labels=c("Reached Vote", "Became Law", "Prop Yeas", "Prop Minority Yeas", "Unity vote"),
          title='Estimates and standard errors for models showing the relationship between veto override threshold and our dependent variables. These results are displayed graphically in Figure 3 of the main text.',
          covariate.labels = c("Veto Override Threshold", "Divided Government", "Line-Item Veto", "Majority Sponsor", 
                               "Number of Cosponsors", "Majority Seat Share", "Party Difference", "Legislative Professionalism",
                               "Committee Sponsor", "Majority Unity", "Governor's Party can Override",
                               2011:2014, "Veto Override Threshold * Divided Government"))

# Term difference
term_mods <- list(mods_vote_divided_bicam_int[[2]], mods_law_divided_bicam_int[[2]], mods_yea_divided_bicam_int[[2]], 
                  mods_yea_min_divided_bicam_int[[2]], mods_unity_divided_bicam_int[[2]])
stargazer(term_mods, keep.stat = 'n', font.size='tiny', column.sep.width='1pt', digits=2, 
          dep.var.caption='', dep.var.labels=c("Reached Vote", "Became Law", "Prop Yeas", "Prop Minority Yeas", "Unity vote"), 
          header=FALSE, label='t.term.bicam', model.names=FALSE, model.numbers=FALSE,
          title='Estimates and standard errors for models showing the relationship between term difference and our dependent variables. These results are displayed graphically in Figure 5 of the main text.',
          covariate.labels = c("Term Difference", 'Divided Government', "Veto Override Threshold", "Line-Item Veto", "Majority Sponsor", 
                               "Number of Cosponsors", "Majority Seat Share", "Party Difference", "Legislative Professionalism", 
                               "Committee Sponsor", "Majority Unity", "Governor's Party can Override", 
                               2011:2014, "Term Difference * Divided Government"))

# Seat ratio
seat_mods <- list(mods_vote_divided_bicam_int[[1]], mods_law_divided_bicam_int[[1]], mods_yea_divided_bicam_int[[1]], 
                  mods_yea_min_divided_bicam_int[[1]], mods_unity_divided_bicam_int[[1]])
stargazer(seat_mods, keep.stat = 'n', font.size='tiny', column.sep.width='1pt', digits=2, 
          dep.var.caption='', dep.var.labels=c("Reached Vote", "Became Law", "Prop Yeas", "Prop Minority Yeas", "Unity vote"), 
          header=FALSE, label='t.seat.bicam', model.names=FALSE, model.numbers=FALSE,
          title='Estimates and standard errors for models showing the relationship between seat ratio and our dependent variables. These results are displayed graphically in Figure 6 of the main text.',
          covariate.labels = c("Seat Ratio", 'Divided Government', "Veto Override Threshold", "Line-Item Veto", "Majority Sponsor", 
                               "Number of Cosponsors", "Majority Seat Share", "Party Difference", "Legislative Professionalism", 
                               "Committee Sponsor", "Majority Unity", "Governor's Party can Override", 
                               2011:2014, "Term Difference * Divided Government"))



###############################################################################################################################################################
### Restricted veto override threshold regression table (Table 7)

# Dataframes
thresh <- filter(votes, (yea_prop > (3/5)) & (yea_prop < (2/3)) & (veto_override_threshold!=0.5) ) %>%
  mutate(low_threshold = ifelse(veto_override_threshold==0.6, '3/5', '2/3') %>% as.factor() %>% relevel(ref='3/5')) %>%
  dplyr::select(c(bill_id, law_sles, low_threshold, divided3, maj_sponsor, num_cs, maj_share, pty_diff, mds1, chamber, committee_sponsor, maj_unity, gov_pty_override, lineitemveto, year, state)) %>% 
  distinct()
thresh2 <- filter(votes, (yea_prop > (1/2)) & (yea_prop < (3/5)) & (veto_override_threshold!=0.67) ) %>%
  mutate(low_threshold = ifelse(veto_override_threshold==0.5, '1/2', '3/5') %>% as.factor() %>% relevel(ref='1/2')) %>%
  dplyr::select(c(bill_id, law_sles, low_threshold, divided3, maj_sponsor, num_cs, maj_share, pty_diff, mds1, chamber, committee_sponsor, maj_unity, gov_pty_override, lineitemveto, year, state)) %>% 
  distinct()
thresh3 <- filter(votes, (yea_prop > (1/2)) & (yea_prop < (2/3)) & (veto_override_threshold!=0.6) ) %>%
  mutate(low_threshold = ifelse(veto_override_threshold==0.5, '1/2', '2/3') %>% as.factor() %>% relevel(ref='1/2')) %>%
  dplyr::select(c(bill_id, law_sles, low_threshold, divided3, maj_sponsor, num_cs, maj_share, pty_diff, mds1, chamber, committee_sponsor, maj_unity, gov_pty_override, lineitemveto, year, state)) %>% 
  distinct()

# Models
mods_thresh <- glm(law_sles ~ low_threshold * divided3 + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + chamber + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(thresh, !state %in% high_missingness), family='binomial')
mods_thresh2 <- glm(law_sles ~ low_threshold * divided3 + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + chamber + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(thresh2, !state %in% high_missingness), family='binomial')
mods_thresh3 <- glm(law_sles ~ low_threshold * divided3 + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + chamber + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(thresh3, !state %in% high_missingness), family='binomial')

# Regression table
thresh_mods <- list(mods_thresh, mods_thresh2, mods_thresh3)
stargazer(thresh_mods, keep.stat='n', column.sep.width = '1pt', digits=2, dep.var.caption = '', dep.var.labels = 'Became Law',
          header=FALSE, label='t.thresh.focus', model.names=FALSE, model.numbers = FALSE, title='', no.space=TRUE, font.size='footnotesize',
          covariate.labels = c("Lower Threshold", "Lower Threshold", "Divided Government", "Majority Sponsor", "Number of Cosponsors", 
                               "Majority Seat Share", "Party Difference", "Legislative Professionalism", 
                               "Senate", "Committee Sponsor", "Majority Unity", "Governor's Party can Override", 2011:2014,
                               "Lower Threshold * Divided Government", "Lower Threshold * Divided Government"))



###############################################################################################################################################################
### Nebraska analysis regression table (Table 8)

# Run models
mods_unicam_vote <- vector(mode='list', length=2)
mods_unicam_vote[[1]] <- glm(voted_sles ~ bicameral + num_cs + as.factor(year), data=bicam)
mods_unicam_vote[[2]] <- glm(voted_sles ~ bicameral + num_cs + as.factor(year), data=bicam2)

mods_unicam_law <- vector(mode='list', length=2)
mods_unicam_law[[1]] <- glm(law_sles ~ bicameral + num_cs + as.factor(year), data=bicam)
mods_unicam_law[[2]] <- glm(law_sles ~ bicameral + num_cs + as.factor(year), data=bicam2)

# Make table
bicam_mods <- list(mods_unicam_vote, mods_unicam_law)
stargazer(bicam_mods, keep.stat='n', column.sep.width = '1pt', digits=2, dep.var.caption = '',
          dep.var.labels = c("Reached Vote", "Became Law"),
          header=FALSE, label='t.bicam.itself', model.names=FALSE, model.numbers = FALSE,
          title='',
          covariate.labels = c("Bicameralism (NE)", "Number of Cosponsors", "2011", "2012", "2013", "2014"))



########################################################################################################################################################################
### Including high-missingness states (Figures 2-4)

### Proportion of members voting yea

# Run models
mods_yea_no_miss <- vector(mode='list', length=5)
mods_yea_no_miss[[1]] <- lm(yea_prop ~ veto_override_threshold*divided + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)
mods_yea_no_miss[[2]] <- lm(yea_prop ~ lineitemveto_plot*divided + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)
mods_yea_no_miss[[3]] <- lm(yea_prop ~ seat_ratio*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)
mods_yea_no_miss[[4]] <- lm(yea_prop ~ term_diff2*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)
mods_yea_no_miss[[5]] <- lm(yea_prop ~ pty_diff*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)

# Predicted effects 
plot_yea_no_miss_veto <- plot_model(mods_yea_no_miss[[1]], type='pred', terms=c('veto_override_threshold', 'divided'), colors='bw') +
  labs(x='Veto Override Threshold', y='Proportion of Members Voting Yea', title='', linetype='') + theme_bw()
plot_yea_no_miss_line <- plot_model(mods_yea_no_miss[[2]], type='pred', terms=c('lineitemveto_plot', 'divided'), colors='bw') +
  labs(x='Line Item Veto', y='Proportion of Members Voting Yea', title='', shape='') + theme_bw()
plot_yea_no_miss_seat <- plot_model(mods_yea_no_miss[[3]], type='pred', terms=c('seat_ratio', 'divided'), colors='bw') +
  labs(x='Seat Ratio', y='Proportion of Members Voting Yea', title='', linetype='') + theme_bw()
plot_yea_no_miss_term <- plot_model(mods_yea_no_miss[[4]], type='pred', terms=c('term_diff2', 'divided'), colors='bw') +
  labs(x='Difference in Term Length between Chambers', y='Proportion of Members Voting Yea', title='', shape='') + theme_bw()
plot_yea_no_miss_pty_diff <- plot_model(mods_yea_no_miss[[5]], type='pred', terms=c('pty_diff', 'divided'), colors='bw') +
  labs(x='Absolute Difference between Party Medians', y='Proportion of Members Voting Yea', title='', linetype='') + theme_bw()

# Put together, save
plot_yea_divided <- (plot_yea_no_miss_veto|plot_yea_no_miss_line)/(plot_yea_no_miss_seat + plot_yea_no_miss_term + plot_yea_no_miss_pty_diff) + 
  plot_annotation(title='Proportion of Members Voting yea', theme = theme(plot.title = element_text(hjust = 0.5)))
ggsave('results/yea_no_miss.pdf', plot=plot_yea_divided, width=10, height=10)


### Proportion of minority members voting yea

# Run models
mods_yea_min_no_miss <- vector(mode='list', length=5)
mods_yea_min_no_miss[[1]] <- lm(min_yea_prop ~ veto_override_threshold*divided + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)
mods_yea_min_no_miss[[2]] <- lm(min_yea_prop ~ lineitemveto_plot*divided + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)
mods_yea_min_no_miss[[3]] <- lm(min_yea_prop ~ seat_ratio*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)
mods_yea_min_no_miss[[4]] <- lm(min_yea_prop ~ term_diff2*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)
mods_yea_min_no_miss[[5]] <- lm(min_yea_prop ~ pty_diff*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, data=votes)

# Predicted effects 
plot_yea_min_no_miss_veto <- plot_model(mods_yea_min_no_miss[[1]], type='pred', terms=c('veto_override_threshold', 'divided'), colors='bw') +
  labs(x='Veto Override Threshold', y='Proportion of Minority Members Voting Yea', title='', linetype='') + theme_bw()
plot_yea_min_no_miss_line <- plot_model(mods_yea_min_no_miss[[2]], type='pred', terms=c('lineitemveto_plot', 'divided'), colors='bw') +
  labs(x='Line Item Veto', y='Proportion of Minority Members Voting Yea', title='', shape='') + theme_bw()
plot_yea_min_no_miss_seat <- plot_model(mods_yea_min_no_miss[[3]], type='pred', terms=c('seat_ratio', 'divided'), colors='bw') +
  labs(x='Seat Ratio', y='Proportion of Minority Members Voting Yea', title='', linetype='') + theme_bw()
plot_yea_min_no_miss_term <- plot_model(mods_yea_min_no_miss[[4]], type='pred', terms=c('term_diff2', 'divided'), colors='bw') +
  labs(x='Difference in Term Length between Chambers', y='Proportion of Minority Members Voting Yea', title='', shape='') + theme_bw()
plot_yea_min_no_miss_pty_diff <- plot_model(mods_yea_min_no_miss[[5]], type='pred', terms=c('pty_diff', 'divided'), colors='bw') +
  labs(x='Absolute Difference between Party Medians', y='Proportion of Minority Members Voting Yea', title='', linetype='') + theme_bw()

# Put together, save
plot_yea_min_divided <- (plot_yea_min_no_miss_veto|plot_yea_min_no_miss_line)/(plot_yea_min_no_miss_seat + plot_yea_min_no_miss_term + plot_yea_min_no_miss_pty_diff) + 
  plot_annotation(title='Proportion of Minority Members Voting Yea', theme = theme(plot.title = element_text(hjust = 0.5)))
ggsave('results/yea_min_min_no_miss.pdf', plot=plot_yea_min_divided, width=10, height=10)


### Pr(Unity Vote)

# Run models
mods_unity_no_miss <- vector(mode='list', length=5)
mods_unity_no_miss[[1]] <- glm(unity_vote ~ veto_override_threshold*divided + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, family='binomial', data=votes)
mods_unity_no_miss[[2]] <- glm(unity_vote ~ lineitemveto_plot*divided + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, family='binomial', data=votes)
mods_unity_no_miss[[3]] <- glm(unity_vote ~ seat_ratio*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, family='binomial', data=votes)
mods_unity_no_miss[[4]] <- glm(unity_vote ~ term_diff2*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, family='binomial', data=votes)
mods_unity_no_miss[[5]] <- glm(unity_vote ~ pty_diff*divided + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + maj_pass_override + year + committee_sponsor, family='binomial', data=votes)

# Predicted effects 
plot_unity_no_miss_veto <- plot_model(mods_unity_no_miss[[1]], type='pred', terms=c('veto_override_threshold', 'divided'), colors='bw') +
  labs(x='Veto Override Threshold', y='Probability of Party Unity Vote', title='', linetype='') + theme_bw()
plot_unity_no_miss_line <- plot_model(mods_unity_no_miss[[2]], type='pred', terms=c('lineitemveto_plot', 'divided'), colors='bw') +
  labs(x='Line Item Veto', y='Probability of Party Unity Vote', title='', shape='') + theme_bw()
plot_unity_no_miss_seat <- plot_model(mods_unity_no_miss[[3]], type='pred', terms=c('seat_ratio', 'divided'), colors='bw') +
  labs(x='Seat Ratio', y='Probability of Party Unity Vote', title='', linetype='') + theme_bw()
plot_unity_no_miss_term <- plot_model(mods_unity_no_miss[[4]], type='pred', terms=c('term_diff2', 'divided'), colors='bw') +
  labs(x='Difference in Term Length between Chambers', y='Probability of Party Unity Vote', title='', shape='') + theme_bw()
plot_unity_no_miss_pty_diff <- plot_model(mods_unity_no_miss[[5]], type='pred', terms=c('pty_diff', 'divided'), colors='bw') +
  labs(x='Absolute Difference between Party Medians', y='Probability of PartyUnity Vote', title='', linetype='') + theme_bw()

# Put together, save
plot_unity_divided <- (plot_unity_no_miss_veto|plot_unity_no_miss_line)/(plot_unity_no_miss_seat + plot_unity_no_miss_term + plot_unity_no_miss_pty_diff) + 
  plot_annotation(title='Probability of Party Unity Vote', theme = theme(plot.title = element_text(hjust = 0.5)))
ggsave('results/unity_no_miss.pdf', plot=plot_unity_divided, width=10, height=10)



########################################################################################################################################################################
### NC CODE GOES HERE


# PUT THAT CODE HERE


########################################################################################################################################################################
### Removing minority-sponsored bills (Figures 7-10)

### Bicameralism models

# Seat ratio
mods_seat_ratio <- vector(mode='list', length=5)
mods_seat_ratio[[1]] <- glm(voted_sles ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_seat_ratio[[2]] <- glm(law_sles ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_seat_ratio[[3]] <- lm(yea_prop ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_seat_ratio[[4]] <- lm(min_yea_prop ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_seat_ratio[[5]] <- glm(unity_vote ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))

# Term limits
mods_term <- vector(mode='list', length=5)
mods_term[[1]] <- glm(voted_sles ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_term[[2]] <- glm(law_sles ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_term[[3]] <- lm(yea_prop ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_term[[4]] <- lm(min_yea_prop ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_term[[5]] <- glm(unity_vote ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))


### Bicameralism plots

# Set up labels
ylabs <- c("Probability of Reaching a Vote", "Probability of Becoming Law", "Proportion of Members Voting Yea",
           "Proportion of Minority Members Voting Yea", "Probability of Party Unity Vote")

# Empty lists to hold plots
plots_seat <- vector(mode='list', length=d(mods_seat_ratio)+1)
plots_term <- vector(mode='list', length=d(mods_seat_ratio)+1)

# Loop through, make plots
for(ii in 1:d(mods_seat_ratio)){
  plots_seat[[ii]] <- plot_basic(mods_seat_ratio[[ii]], c("seat_ratio", "divided3"), "Seat Ratio", ylabs[ii])
  plots_term[[ii]] <- plot_basic(mods_term[[ii]], c("term_diff2", "divided3"), "Term Difference", ylabs[ii])
}

# Get legends
seat_leg <- plot_model(mods_seat_ratio[[1]], type='pred', line.size=2,
                       terms=c('seat_ratio', 'divided3'), colors='bw') +
  labs(linetype='Legend') + theme_bw() +
  theme(legend.text=element_text(size=20), legend.title=element_text(size=20)) +
  guides(linetype = guide_legend(override.aes = list(size = 20)))
plots_seat[[6]] <- cowplot::get_legend(seat_leg) %>%
  as_ggplot() + theme(legend.key.size=unit(10,'cm'))

term_leg <- plot_model(mods_term[[1]], type='pred', line.size=2,
                       terms=c('term_diff2', 'divided3'), colors='bw') +
  labs(shape='Legend') + theme_bw() +
  theme(legend.text=element_text(size=20), legend.title=element_text(size=20)) +
  guides(shape = guide_legend(override.aes = list(size = 15)))
plots_term[[6]] <- cowplot::get_legend(term_leg) %>%
  as_ggplot() + theme(legend.key.size=unit(10,'cm'))

# Plot together
seat_plots_combined <- wrap_plots(plots_seat) + plot_layout(ncol=2)
term_plots_combined <- wrap_plots(plots_term) + plot_layout(ncol=2)

# Save
ggsave('results/seat_divided_maj_only_plot.pdf', width=12, height=15, plot=seat_plots_combined)
ggsave('results/term_divided_maj_only_plotplot.pdf', width=12, height=15, plot=term_plots_combined)


### Interchamber models

# Seat ratio
mods_veto <- vector(mode='list', length=5)
mods_veto[[1]] <- glm(voted_sles ~ veto_override_threshold*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_veto[[2]] <- glm(law_sles ~ veto_override_threshold*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_veto[[3]] <- lm(yea_prop ~ veto_override_threshold*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_veto[[4]] <- lm(min_yea_prop ~ veto_override_threshold*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_veto[[5]] <- glm(unity_vote ~ veto_override_threshold*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))

# Term limits
mods_line <- vector(mode='list', length=5)
mods_line[[1]] <- glm(voted_sles ~ lineitemveto_plot*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_line[[2]] <- glm(law_sles ~ lineitemveto_plot*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_line[[3]] <- lm(yea_prop ~ lineitemveto_plot*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_line[[4]] <- lm(min_yea_prop ~ lineitemveto_plot*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))
mods_line[[5]] <- glm(unity_vote ~ lineitemveto_plot*divided3 + veto_override_threshold + lineitemveto + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness & maj_sponsor=='Majority'))


### Interchamber plots

# Set up labels
ylabs <- c("Probability of Reaching a Vote", "Probability of Becoming Law", "Proportion of Members Voting Yea",
           "Proportion of Minority Members Voting Yea", "Probability of Party Unity Vote")

# Empty lists to hold plots
plots_veto <- vector(mode='list', length=d(mods_veto)+1)
plots_line <- vector(mode='list', length=d(mods_veto)+1)

# Loop through, make plots
for(ii in 1:d(mods_veto)){
  plots_veto[[ii]] <- plot_basic(mods_veto[[ii]], c("veto_override_threshold", "divided3"), "Veto Override Threshold", ylabs[ii])
  plots_line[[ii]] <- plot_basic(mods_line[[ii]], c("lineitemveto_plot", "divided3"), "Line Item Veto", ylabs[ii])
}

# Get legends
veto_leg <- plot_model(mods_veto[[1]], type='pred', line.size=2,
                       terms=c('veto_override_threshold', 'divided3'), colors='bw') +
  labs(linetype='Legend') + theme_bw() +
  theme(legend.text=element_text(size=20), legend.title=element_text(size=20)) +
  guides(linetype = guide_legend(override.aes = list(size = 20)))
plots_veto[[6]] <- cowplot::get_legend(veto_leg) %>%
  as_ggplot() + theme(legend.key.size=unit(10,'cm'))

line_leg <- plot_model(mods_line[[1]], type='pred', line.size=2,
                       terms=c('lineitemveto_plot', 'divided3'), colors='bw') +
  labs(shape='Legend') + theme_bw() +
  theme(legend.text=element_text(size=20), legend.title=element_text(size=20)) +
  guides(shape = guide_legend(override.aes = list(size = 15)))
plots_line[[6]] <- cowplot::get_legend(line_leg) %>%
  as_ggplot() + theme(legend.key.size=unit(10,'cm'))

# Plot together
veto_plots_combined <- wrap_plots(plots_veto) + plot_layout(ncol=2)
line_plots_combined <- wrap_plots(plots_line) + plot_layout(ncol=2)

# Save
ggsave('results/veto_override_divided_maj_only_plot.pdf', width=12, height=15, plot=veto_plots_combined)
ggsave('results/line_item_divided_maj_only_plot.pdf', width=12, height=15, plot=line_plots_combined)



########################################################################################################################################################################
### Only under divided government (Figures 11-14)

### Bicameralism models

# Seat ratio
mods_seat_ratio <- vector(mode='list', length=5)
mods_seat_ratio[[1]] <- glm(voted_sles ~ seat_ratio + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & divided3=='Divided'))
mods_seat_ratio[[2]] <- glm(law_sles ~ seat_ratio + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & divided3=='Divided'))
mods_seat_ratio[[3]] <- lm(yea_prop ~ seat_ratio + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & divided3=='Divided'))
mods_seat_ratio[[4]] <- lm(min_yea_prop ~ seat_ratio + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & divided3=='Divided'))
mods_seat_ratio[[5]] <- glm(unity_vote ~ seat_ratio + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness & divided3=='Divided'))

# Term limits
mods_term <- vector(mode='list', length=5)
mods_term[[1]] <- glm(voted_sles ~ term_diff2 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & divided3=='Divided'))
mods_term[[2]] <- glm(law_sles ~ term_diff2 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & divided3=='Divided'))
mods_term[[3]] <- lm(yea_prop ~ term_diff2 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & divided3=='Divided'))
mods_term[[4]] <- lm(min_yea_prop ~ term_diff2 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & divided3=='Divided'))
mods_term[[5]] <- glm(unity_vote ~ term_diff2 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness & divided3=='Divided'))


### Bicameralism plots

# Function 
plot_basic <- function(model, terms, xlab, ylab){
  plot_model(model, type='pred', line.size=2, dot.size=5,
             terms=terms, colors='bw') +
    labs(x=xlab, y=ylab, title='', linetype='') + 
    theme_bw() +
    theme(text=element_text(size=16),
          axis.text=element_text(size=16),
          legend.position="none")
}

# Set up labels
ylabs <- c("Probability of Reaching a Vote", "Probability of Becoming Law", "Proportion of Members Voting Yea",
           "Proportion of Minority Members Voting Yea", "Probability of Party Unity Vote")

# Empty lists to hold plots
plots_seat <- vector(mode='list', length=d(mods_seat_ratio))
plots_term <- vector(mode='list', length=d(mods_seat_ratio))

# Loop through, make plots
for(ii in 1:d(mods_seat_ratio)){
  plots_seat[[ii]] <- plot_basic(mods_seat_ratio[[ii]], c("seat_ratio"), "Seat Ratio", ylabs[ii])
  plots_term[[ii]] <- plot_basic(mods_term[[ii]], c("term_diff2"), "Term Difference", ylabs[ii])
}

# Plot together
seat_plots_combined <- wrap_plots(plots_seat) + plot_layout(ncol=2)
term_plots_combined <- wrap_plots(plots_term) + plot_layout(ncol=2)

# Save
ggsave('results/seat_divided_only_plot.pdf', width=12, height=15, plot=seat_plots_combined)
ggsave('results/term_divided_only_plot.pdf', width=12, height=15, plot=term_plots_combined)


### Interchamber models

# Seat ratio
mods_veto <- vector(mode='list', length=5)
mods_veto[[1]] <- glm(voted_sles ~ veto_override_threshold + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & divided3=='Divided'))
mods_veto[[2]] <- glm(law_sles ~ veto_override_threshold + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & divided3=='Divided'))
mods_veto[[3]] <- lm(yea_prop ~ veto_override_threshold + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & divided3=='Divided'))
mods_veto[[4]] <- lm(min_yea_prop ~ veto_override_threshold + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & divided3=='Divided'))
mods_veto[[5]] <- glm(unity_vote ~ veto_override_threshold + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness & divided3=='Divided'))

# Term limits
mods_line <- vector(mode='list', length=5)
mods_line[[1]] <- glm(voted_sles ~ lineitemveto_plot + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & divided3=='Divided'))
mods_line[[2]] <- glm(law_sles ~ lineitemveto_plot + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness & divided3=='Divided'))
mods_line[[3]] <- lm(yea_prop ~ lineitemveto_plot + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & divided3=='Divided'))
mods_line[[4]] <- lm(min_yea_prop ~ lineitemveto_plot + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness & divided3=='Divided'))
mods_line[[5]] <- glm(unity_vote ~ lineitemveto_plot + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness & divided3=='Divided'))


### Interchamber plots

# Set up labels
ylabs <- c("Probability of Reaching a Vote", "Probability of Becoming Law", "Proportion of Members Voting Yea",
           "Proportion of Minority Members Voting Yea", "Probability of Party Unity Vote")

# Empty lists to hold plots
plots_veto <- vector(mode='list', length=d(mods_veto))
plots_line <- vector(mode='list', length=d(mods_veto))

# Loop through, make plots
for(ii in 1:d(mods_veto)){
  plots_veto[[ii]] <- plot_basic(mods_veto[[ii]], c("veto_override_threshold"), "Veto Override Threshold", ylabs[ii])
  plots_line[[ii]] <- plot_basic(mods_line[[ii]], c("lineitemveto_plot"), "Line Item Veto", ylabs[ii])
}

# Plot together
veto_plots_combined <- wrap_plots(plots_veto) + plot_layout(ncol=2)
line_plots_combined <- wrap_plots(plots_line) + plot_layout(ncol=2)

# Save
ggsave('results/veto_override_divided_only_plot.pdf', width=12, height=15, plot=veto_plots_combined)
ggsave('results/line_item_divided_only_plot.pdf', width=12, height=15, plot=line_plots_combined)



##########################################################################################################################################################
### Comparing line-item veto effects among states with similar institutional arrangements (Figure 15)

### Format data: get states that otherwise match the institutional setup for at least one of the no lineitemveto states (NOT including seat ratio)

# Get state institutions data
st_insts <- bill %>%
  filter(!is.na(lineitemveto) & !is.na(veto_override_threshold) & !is.na(seat_ratio) & !is.na(term_diff) & state!='NE') %>%
  group_by(state) %>%
  dplyr::reframe(
    lineitemveto = unique(lineitemveto),
    veto_override_threshold = unique(veto_override_threshold),
    seat_ratio = unique(seat_ratio),
    term_diff = unique(term_diff)
  ) %>%
  ungroup() %>%
  filter(!(state=='WY' & veto_override_threshold==0.5)) %>% # Wyoming has different thresholds by chamber, so functionally there is a 0.67 threshold for passage, so we'll use that
  filter(!(state=='NC' & lineitemveto==1)) # Typo in NC data for 2010

# Pull out the profile for non-lineitemveto states
no_liv <- filter(st_insts, lineitemveto==0)

# Go through each state, check if it matches any of those in no_liv
state_matches <- c()
for(ii in 1:nrow(st_insts)){
  if(!st_insts$state[ii] %in% no_liv$state){
    match_veto <- which(no_liv$veto_override_threshold == st_insts$veto_override_threshold[ii])
    match_term <- which(no_liv$term_diff == st_insts$term_diff[ii])
    # match_seat <- (no_liv$seat_ratio-st_insts$seat_ratio[ii]) %>%
    #   abs() %>%
    #   `<`(0.1) %>%
    #   which()
    if(d(intersect(match_veto, match_term)) > 0){
      res <- st_insts$state[ii]
      state_matches <- c(state_matches, res)
    }
  }
}

# Make dataframes
match_bill <- filter(bill, state!='NE' & (lineitemveto==0 | state %in% state_matches))
match_vote <- filter(votes, state!='NE' & (lineitemveto==0 | state %in% state_matches))


### Format data: get states that match on veto override threshold, term length difference, *and* seat ratio (within 0.1); only need to compare to NV and VT because the others don't have a close match for seat ratio

# List of states matching NV or VT
# filter(st_insts, state %in% c("NV", "VT"))
nv_match <- filter(st_insts, veto_override_threshold==0.67 & abs(seat_ratio-2) < 0.1 & term_diff==2 & lineitemveto==1) %>%
  pull(state)
vt_match <- filter(st_insts, veto_override_threshold==0.67 & abs(seat_ratio-1.25) < 0.1 & term_diff==0 & lineitemveto==1) %>%
  pull(state)

# Combine
line_match <- c(nv_match, vt_match, "NV", "VT") %>% unique()

# Make dataframes
match_bill_sml <- filter(bill, state %in% line_match)
match_vote_sml <- filter(votes, state %in% line_match)


### Models without considering a match in seat_ratio

# Run models
mod_vote <- glm(voted_sles ~ lineitemveto_plot*divided3 + veto_override_threshold + term_diff + seat_ratio + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=match_bill)
mod_law <- glm(law_sles ~ lineitemveto_plot*divided3 + veto_override_threshold + term_diff + seat_ratio + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=match_bill)
mod_yea <- lm(yea_prop ~ lineitemveto_plot*divided3 + veto_override_threshold + term_diff + seat_ratio + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=match_vote)
mod_yea_min <- lm(min_yea_prop ~ lineitemveto_plot*divided3 + veto_override_threshold + term_diff + seat_ratio + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=match_vote)
mod_unity <- lm(unity_vote ~ lineitemveto_plot*divided3 + veto_override_threshold + term_diff + seat_ratio + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=match_vote)

# Substantive effect sizes
sub_eff_vote <- sub_eff(model=mod_vote, terms=c("lineitemveto_plot", "divided3"), dv="voted_sles", data=match_bill)
sub_eff_law <- sub_eff(model=mod_law, terms=c("lineitemveto_plot", "divided3"), dv="law_sles", data=match_bill)
sub_eff_yea <- sub_eff(model=mod_yea, terms=c("lineitemveto_plot", "divided3"), dv="yea_prop", data=match_vote)
sub_eff_yea_min <- sub_eff(model=mod_yea_min, terms=c("lineitemveto_plot", "divided3"), dv="min_yea_prop", data=match_vote)
sub_eff_unity <- sub_eff(model=mod_unity, terms=c("lineitemveto_plot", "divided3"), dv="unity_vote", data=match_vote)

# Plots
plot_vote <- plot_eff(model=mod_vote, terms=c("lineitemveto_plot", "divided3"), x="Line-Item Veto",
                      y="Probability of Reaching a Vote", text_df=sub_eff_vote, nudge_x=-0.15)
plot_law <- plot_eff(model=mod_law, terms=c("lineitemveto_plot", "divided3"), x="Line-Item Veto",
                     y="Probability of Becoming Law", text_df=sub_eff_law, nudge_x=-0.15)
plot_yea <- plot_eff(model=mod_yea, terms=c("lineitemveto_plot", "divided3"), x="Line-Item Veto",
                     y="Proportion of Members Voting Yea", text_df=sub_eff_yea, nudge_x=-0.15)
plot_yea_min <- plot_eff(model=mod_yea_min, terms=c("lineitemveto_plot", "divided3"), x="Line-Item Veto",
                         y="Proportion of Minority Members Voting Yea", text_df=sub_eff_yea_min, nudge_x=-0.15)
plot_unity <- plot_eff(model=mod_unity, terms=c("lineitemveto_plot", "divided3"), x="Line-Item Veto",
                       y="Probability of Party Unity Vote", text_df=sub_eff_unity, nudge_x=-0.15)

# Get legend
line_leg <- plot_model(mod_vote, type='pred', dot.size=5,
                       terms=c('lineitemveto_plot', 'divided3'), colors='bw') +
  labs(shape='Legend') + theme_bw() +
  theme(legend.text=element_text(size=20), legend.title=element_text(size=20)) +
  guides(shape = guide_legend(override.aes = list(size = 15))) 
line_leg_plot <- cowplot::get_legend(line_leg) %>%
  as_ggplot() + 
  theme(legend.key.size=unit(10,'cm'))

# Combine
line_plots <- plot_vote + plot_law + plot_yea + plot_yea_min + plot_unity + line_leg_plot + plot_layout(ncol=2)

# Save
ggsave("results/line_item_match_plot.pdf", width=12, height=15, plot=line_plots)
